import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile
from IPython.display import Audio, display
import warnings
warnings.filterwarnings("ignore")
#Fenêtre de hamming
def hamming(T): return 0.54-0.46*np.cos(2*np.pi*np.arange(T)/(T-1))
#Calcul du spectre amplitude et phase du signal data, taille de la fenetre=T, pas=p
#mettre pre=True dans l'appel pour activer la préaccentuation
#mettre ham=True pour activer le fenêtrage par hamming
#Tfft=taille de la fft, si > T le zéro padding sera activé(voir cours acoustique)
def spectrogram(data,fs,T=512,p=32,pre=False,ham=True,Tfft=None,norm=True):
if Tfft is None: Tfft=T #si la taille de la fft n'est pas spécifiée prendre Tfft=T
if norm: data=(data-np.mean(data))/np.std(data) #normaliser le signal sur son écart type
if pre: data[1:]-0.97*data[:-1] #préaccentuation
s=[data[i:i+T] for i in range(0,len(data)-T,p)] # fenêtrage
if ham : s=s*hamming(T) # multiplication par hamming
s=np.fft.fft(s,Tfft) #Transformée de Fourier
s=s[:,:int(Tfft/2)] #couper le spectre en 2 pour éliminer l'effet mirroir
#retourner les spectres d'amplitude et de phase
return {"ampl": np.abs(s), "phase": np.angle(s), "fs":fs, "duree":len(data)/fs, "T": T, "p": p, "Tfft": Tfft}
#Cette fonction affiche un spectrogramme d'amplitude
def afficher(spec):
fig, ax1 = plt.subplots()
ax2 = ax1.twinx()
ax3 = ax1.twiny()
freq=np.linspace(0, spec["fs"]/2, int(spec["Tfft"]/2)) #labels de l'axe des fréquences
temps=np.linspace(0, spec["duree"], len(spec["ampl"])) #labels de l'axe du temps
ax1.pcolormesh(temps, freq, np.log(spec["ampl"]+1).T, cmap='gray_r') #afficher le spectre avec correction des axes
ax1.set_ylabel('Fréquence (Hz)')
ax1.set_xlabel('Temps (sec)')
ax2.set_yticks(np.round(ax2.get_yticks()*int(spec["Tfft"]/2)))
ax2.set_ylabel('Indices des colonnes dans la matrice spec')
ax3.set_xticks(np.round(ax3.get_xticks()*len(spec["ampl"])))
ax3.set_xlabel('Indices des lignes dans la matrice spec')
plt.show()
#Retour à partir de l'amplitude et de la phase au signal temporel
#Ne pas utiliser la préaccentuation si vous voulez utiliser cette fonction
def spec2wav(spec):
p=spec["p"]; T=spec["T"]
ne=(len(spec["ampl"])-1)*p+T #estimation le nombre d'échantillons du signal
signal=np.zeros(ne) #initialiser le signal par des zéros
trams=np.zeros(ne) #initialiser le nombre de trames par des zéros
temp=spec["ampl"]*np.exp(1j*spec["phase"]) #recombiner l'amplitude avec la phase (nombre complexe)
temp=np.fft.ifft(temp, spec["Tfft"]) #retour au domaine temporel par une fft inverse
temp=np.real(temp) #ne garder que la partie réelle
for i in range(len(temp)): #fenêtrage inverse
signal[i*p:i*p+T]+=temp[i,:T]
trams[i*p:i*p+T]+=1
return signal/trams, spec["fs"] #retourner le signal reconstitué et la fréquence d'échantillonnage
dataset_path="Desktop\\TP dataset\\dataset\\" #le chemin vers le dataset
filename="7_george_0.wav" #le nom du fichier à charger
print('Chargement de :', filename)#chargement d'un fichier audio wav
fs, data = wavfile.read(dataset_path+filename) #fs: fréquence d'échantillonnage
print('Freq échantillonnage:',fs,'Hz, durée:',len(data)/fs,'sec')
display(Audio(data,rate=fs))
plt.plot(np.arange(len(data))/fs,data); plt.show() #affichage du signal temporel
spec=spectrogram(data, fs, T=256, Tfft=1025) #calcul du spectrogramme
afficher(spec) #affichage du spectrogramme
Chargement de : 7_george_0.wav Freq échantillonnage: 8000 Hz, durée: 0.641375 sec
newdata, fs= spec2wav(spec)
display(Audio(newdata,rate=fs))
spec=spectrogram(data, fs, T=256, Tfft=1024) #calcul et affichage d'un spectrogramme avec fenêtre de hamming
spec["ampl"][:,128:]=0 #annuler toutes les fréquences > 1000Hz
afficher(spec) #affichage du spectrogramme
newdata, fs= spec2wav(spec) #retour au domaine temporel
display(Audio(newdata,rate=fs)) #écouter le résultat
spec=spectrogram(data, fs, T=256, Tfft=1024) #calcul et affichage d'un spectrogramme avec fenêtre de hamming
spec["ampl"][:,:256]=0 #annuler toutes les fréquences < 2000Hz
afficher(spec) #affichage du spectrogramme
newdata, fs= spec2wav(spec) #retour au domaine temporel
display(Audio(newdata,rate=fs)) #écouter le résultat
spec=spectrogram(data, fs, T=256, Tfft=1024) #calcul et affichage d'un spectrogramme avec fenêtre de hamming
spec["ampl"][:,:128]=0 #annuler toutes les fréquences < 1000Hz
spec["ampl"][:,256:]=0 #annuler toutes les fréquences > 2000Hz
afficher(spec) #affichage du spectrogramme
newdata, fs= spec2wav(spec) #retour au domaine temporel
display(Audio(newdata,rate=fs)) #écouter le résultat
Q1. Essayez les filtres en haut sur les mots two, four, seven, eight et notez les phonèmes qui sont affectés par le filtrage.
Q2. Chargez le fichier mots_bruit.wav et essayez d'y filtrer le bruit avec la méthode de soustraction spectrale.
fs, data = wavfile.read("mots_bruit.wav") #chargement d'un fichier audio wav
print('Signal orginal: Freq échantillonnage:',fs,'Hz, durée:',len(data)/fs,'sec')
display(Audio(data,rate=fs))
Signal orginal: Freq échantillonnage: 8000 Hz, durée: 19.569875 sec
from scipy.interpolate import interp1d
# Etirer le spectre d'amplitude avec le facteur "a" (new f0= f0*a)
def interpolate(spec, a):
n,m=spec.shape
newspec=np.zeros((n,m))
for i in range(n):
f=interp1d(np.arange(m), spec[i,:])
new_x=np.arange(m)/a
new_x[new_x>m-1]=m-1
newspec[i,:]=f(new_x)
return newspec
######## début du programme ####################
fs, data = wavfile.read("fr.wav")
print('Signal orginal: Freq échantillonnage:',fs,'Hz, durée:',len(data)/fs,'sec')
display(Audio(data,rate=fs))
a=0.7
spec=spectrogram(data, fs, p=16, T=512)
spec["ampl"]=interpolate(spec["ampl"], a)
newdata, fs= spec2wav(spec)
print('Signal avec f0 manipulé de %f'%a)
display(Audio(newdata,rate=fs))
a=1.3
spec=spectrogram(data, fs, p=16, T=512)
spec["ampl"]=interpolate(spec["ampl"], a)
newdata, fs= spec2wav(spec)
print('Signal avec f0 manipulé de %f'%a)
display(Audio(newdata,rate=fs))
Signal orginal: Freq échantillonnage: 16000 Hz, durée: 19.93 sec
Signal avec f0 manipulé de 0.700000
Signal avec f0 manipulé de 1.300000
Q1. Essayez la manipulation sur les fichiers fr.wav, ar.wav, en.wav avec différents paramètres "a" et déduire l'impact de la valeur de "a" sur la parole
## Fonctions utiles de conversion
def Mel2Hz(mel): return 700 * (np.power(10,mel/2595)-1)
def Hz2Mel(freq): return 2595 * np.log10(1+freq/700)
def Hz2Ind(freq,fs,Tfft): return (freq*Tfft/fs).astype(int)
#Réalisation de nf filtres sur l'échelle Mel d'une fréquence min à une fréquence max
def FiltresMel(fs, nf, Tfft, fmin, fmax):
Indices=Hz2Ind(Mel2Hz(np.linspace(Hz2Mel(fmin), Hz2Mel(min(fmax,fs/2)), nf+2)),fs,Tfft)
filtres=np.zeros((int(Tfft/2), nf))
for i in range(nf):
f=Indices[i+2]
if (Indices[i]-f)%2==0: f=min(f+1,int(Tfft/2))
if (f-Indices[i])>1: filtres[Indices[i]:f,i]=hamming(f-Indices[i])
else: filtres[Indices[i]:f,i]=1
return filtres
nf=20 #nombre de filtres
T=256
Tfft=1024
fs=8000
#réalisation des filtres auditifs
filtres=FiltresMel(fs, nf=nf, Tfft=Tfft, fmin=500, fmax=3500)
plt.plot(filtres)
#affichage des filtres
plt.xticks(plt.xticks()[0][1:-1],np.round(plt.xticks()[0][1:-1]*fs/Tfft))
plt.xlabel("Fréquence (Hz)")
plt.title("%d Filtres auditifs"%nf)
plt.show()
## Filtrage
filename="7_george_0.wav"
fs, data = wavfile.read(dataset_path+filename)
spec=spectrogram(data, fs, T=T, Tfft=Tfft) #calcul du spectrogramme
afficher(spec); #affichage du spectrogramme
#affichage de la sortie des filtres auditifs
spec_filtre=(spec["ampl"]@filtres)/10
plt.pcolormesh(np.log(spec_filtre+1).T, cmap="gray_r")
plt.title("Spectrogramme filtré par %d filtres auditifs"%nf)
plt.ylabel("Sortie de chaque filtre")
plt.show()
Objectif: Insérer un maximum de zéros en éliminant les sons qui ne sont pas perçus par l'oreille humaine
nf=36 #nombre de filtres
T=512
fs, data = wavfile.read("fr.wav") #chargement d'un fichier audio wav
#Foncrion Compression du signal par les filtres auditifs
def Compression(spec, nf, fmin=200, fmax=4000):
filtres=FiltresMel(spec["fs"], nf, spec["Tfft"], fmin, fmax)
newspec=np.zeros(spec["ampl"].shape)
for i in range(filtres.shape[1]):
ind=np.argmax(spec["ampl"]*(filtres[:,i]>0),axis=1)
for j in range(len(ind)): #ne grader que la valeur max de chaque bande de filtre
newspec[j, ind[j]]=spec["ampl"][j, ind[j]]
return newspec
################# début du programme #################
print('Signal orginal: Freq échantillonnage:',fs,'Hz, durée:',len(data)/fs,'sec')
display(Audio(data,rate=fs))
spec=spectrogram(data, fs, p=int(T/32), T=T) #calcul du spectrogramme
afficher(spec)
plt.show()
newspec=Compression(spec, nf=nf, fmin=0, fmax=8000) #compression
spec["ampl"]=newspec #remplacer le spectre d'amplitude par le spectre compressé
newdata, fs= spec2wav(spec) #retour au domaine temporel
print("Signal compressé avec %d filtres"%nf)
display(Audio(newdata,rate=fs)) #écouter le résultat
afficher(spec)
Signal orginal: Freq échantillonnage: 16000 Hz, durée: 19.93 sec
Signal compressé avec 36 filtres
Sachant que plus le nombre de filtres est petit, plus le son est compressé.
Q1. Essayez plusieurs paramètres nf, T, fmin, fmax, sur les fichier fr.wav, ar.wav, en.wav afin d'obtenir un bon taux de compression avec une bonne qualité de son.